library(rgdal)
library(readr)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(rgdal)
library(tmap)
library(sf)
library(cowplot)
library(lubridate)

shp<- readOGR("data/shapefiles/iraq_districts_06162010.shp")

survmapvars <- read_csv("data/raw/survmapvars.csv")

colnames(survmapvars)[colnames(survmapvars)=="district"] <- "ADM3NAME"

shp <- merge(shp, survmapvars, by="ADM3NAME")

shp@data$survcat <- NA
shp@data$survcat[shp@data$cgp==1] <- 1
shp@data$survcat[shp@data$tgp==1] <- 2
shp@data$survcat[shp@data$bctgp==1] <- 3

shp@data$survcat[shp@data$survcat==1] <- "Control group (before Mosul)"
shp@data$survcat[shp@data$survcat==2] <- "Treatment group (after Mosul)"
shp@data$survcat[shp@data$survcat==3] <- "Treatment group and control group"

shp_data <- fortify(shp, region = "ADM3NAME")
shp_data$id <- factor(shp_data$id)
shp@data$id <- as.factor(shp@data$ADM3NAME)
plot.data <- left_join(shp_data, shp@data, by = "id")
plot.data$survcat <- as.factor(plot.data$survcat)
#final survey respondents map
#add Mosul label
moscd = data.frame(long=c(3993441),
                 lat=c(317004.3),
                 label=c("Mosul"))

smap <- ggplot(plot.data, aes(x = long, y = lat, 
                      group = group, fill = survcat)) +
  geom_polygon(color = "black", size = 0.05) +
  geom_point(data = moscd, aes(x = lat, y = long,
                               shape =  label),
             lwd= 8,inherit.aes = F, col = "black") +
  coord_equal() + 
  theme_tufte(base_family = "Helvetica") +
  scale_shape_manual("", values=c(8)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank(),
        legend.position = "left",
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  theme(legend.position = "left",
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  scale_fill_manual(
    values = c("#BDBDBC", "#5E5E5B", "#000000"),
    na.value = c("#ffffff"),
    name = "Surveyed status",
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      nrow=4,
      byrow=TRUE
    ))
smap

png("plots/survmap.png", width = 465, height = 225, units='mm', res = 300)
smap
dev.off()

# clan ties map

shp <- readOGR("data/shapefiles/iraq_districts_06162010.shp")
clanties <- read.csv("data/raw/clanties.csv")

clanties$clan_ties <- clanties$clan_ties*100
shp <- merge(shp, clanties, by="ADM3NAME")
shp_data <- fortify(shp, region = "ADM3NAME")

shp@data$id <- as.factor(shp@data$ADM3NAME)
plot.data <- left_join(shp_data, shp@data, by = "id")

cmap <- ggplot(plot.data, aes(x = long, y = lat, 
                      group = group, fill =  clan_ties)) +
  geom_polygon(color = "black", size = 0.05) +
  coord_equal() + 
  theme_tufte(base_family = "Helvetica") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank(),
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  guides(fill=guide_legend(title="Tribal ties")) +
  scale_fill_gradient(
    low = "#BDBDBC",
    high = "#000000",
    na.value = "#ffffff"
    )
cmap

png("plots/clansmap.png", width = 465, height = 225, units='mm', res = 300)
cmap
dev.off()

######## ISIS path map

shp<- readOGR("data/shapefiles/iraq_districts_06162010.shp")

survmapvars <- read.csv("data/raw/survmapvars.csv")

colnames(survmapvars)[colnames(survmapvars)=="district"] <- "ADM3NAME"

shp <- merge(shp, survmapvars, by="ADM3NAME")

qtm(shp)
shpsf <- st_as_sf(shp)
qtm(shpsf)

#42.9356,33.6231,44.8912,35.0518
#43.4792,33.2768,45.0777,34.8917

shpsf$tcgp <- NA
shpsf$tcgp[shpsf$cgp==1] <- 0
shpsf$tcgp[shpsf$tgp==1] <- 1
shpsf$tcgp[shpsf$bctgp==1] <- 2

shpsf$tcgpfac <- as.factor(shpsf$tcgp)

shpsf <- st_transform(shpsf, "+proj=longlat +datum=WGS84")

shpsf$tcgp[shpsf$tcgp==0] <- "Control group (before Mosul)"
shpsf$tcgp[shpsf$tcgp==1] <- "Treatment group (after Mosul)"
shpsf$tcgp[shpsf$tcgp==2] <- "Treatment group and control group"

g1 <- ggplot(data = shpsf) +
  geom_sf(aes(fill = tcgp),lwd=.1) +
  geom_rect(xmin = 43, xmax = 45, ymin = 33.4, ymax = 35, 
            fill = NA, colour = "black", size = 1) +
  theme_tufte(base_family = "Helvetica") +
  annotate(geom = "text", x = 43.12, y = 36.6, label = "Mosul",
           color = "black", size = 6) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank(),
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  theme(legend.position = "left",
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  scale_fill_manual(
    values = c("#BDBDBC", "#5E5E5B", "#000000"),
    na.value = c("#ffffff"),
    name = "Surveyed status",
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      nrow=4,
      byrow=TRUE
    ))

g1


g2 <- ggplot(data = shpsf) +
  geom_sf(aes(fill = tcgp),lwd=.1) +
  geom_rect(xmin = 43, xmax = 45, ymin = 33.4, ymax = 35, 
            fill = NA, colour = "black", size = 1) +
  coord_sf(xlim = c(43,45), ylim = c(33.4,35), expand = FALSE) +
  annotate(geom = "text", x = 44.52, y = 34.05, label = "Al-Khalis",
           fontface = "italic", color = "white", size = 5) +
  annotate(geom = "text", x = 44.5, y = 33.55, label = "Ba'quba",
           fontface = "italic", color = "white", size = 4, angle=55) +
  annotate(geom = "text", x = 43.5, y = 34.2, label = "Samarra",
           fontface = "italic", color = "white", size = 6) +
  annotate(geom = "text", x = 43.75, y = 34.8, label = "Tikrit",
           fontface = "italic", color = "white", size = 6) +
  theme_tufte(base_family = "Helvetica") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank(),
        legend.position = "right",
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  theme(legend.position = "none") +
  scale_fill_manual(
    values = c("#BDBDBC", "#5E5E5B", "#000000"),
    na.value = c("#ffffff"))

x1 <- 36.34566
x2 <- 34.0

y1 <- 43.12332
y2 <- 44.0

g1 <- g1 + geom_curve(aes(x = y1, y = x1, xend = y2, yend = x2),
                      curvature = 0.3, col="black", arrow = arrow(length = unit(0.03, "npc"))) +
  annotate(geom = "text", x = 42.3, y = 35.5, label = "ISIS path",
           color = "black", size = 5)

g2 <- g2 +geom_curve(aes(x = y1, y = x1, xend = y2, yend = x2),
                     curvature = 0.3, col="black", arrow = arrow(length = unit(0.03, "npc")))
#add tribes
tribes <- read.csv("data/raw/tribal_mobilization.csv")

g1 <- g1 +  geom_point(data = tribes, aes(x = lat, y = long, 
                                          shape =  location),
                       size=4,
                       stroke = 1.5,
                       col = "black") +
  scale_shape_manual("Location of tribal mobilizations   ", values=c(1,2,3,4,5,6,7,8,9,10,12,13,14,15),
                     guide = guide_legend(
                       direction = "vertical",
                       title.position = "top",
                       legend.position= "left",
                       ncol = 1, size=5
                     ))

g2 <- g2 +  geom_point(data = tribes, aes(x = lat, y = long, 
                                          shape =  location),
                       size=4,
                       stroke = 1.5,
                       col = "black") +
  scale_shape_manual("Location of tribal mobilizations   ", values=c(1,2,3,4,5,6,7,8,9,10,12,13,14,15),
                     guide = guide_legend(
                       direction = "vertical",
                       title.position = "top",
                       legend.position= "left",
                       ncol = 1, size=5
                     ))

png("plots/ispcombined.png", 
    width = 400, height = 200, units='mm', res = 300)
plot_grid(g1, g2, nrow = 1, rel_widths = c(3, 1))
dev.off()

# ---
# Location of survey respondents over time
# ---

pol <- readOGR("data/shapefiles/WGS84/irq_adm_cso_itos_20190606_shp/irq_admbnda_adm1_cso_20190603.shp")
distcoords <- read_csv("data/raw/distesocmatches.csv")
survey_rspdts <- read_csv("data/output/survey_rspdts_times.csv")
survxys <- merge(survey_rspdts, distcoords, by="distat2")

latlon <- distcoords %>%
  select(lat, lon)
y <- latlon$lat
x <- latlon$lon
xy <- cbind(x,y)
xy <- as.data.frame(xy)
pts <- SpatialPointsDataFrame(xy, data = xy)
plot(pol)
points(pts, col="red", pch="+")

#plot all points
coordsts <- survxys %>%
  select(lat, lon, mdyst, timest)
coordsts$mdyst <- as.character(coordsts$mdyst)
coordsts$timest <- as.character(coordsts$timest)
coordsts$dtst <- paste(coordsts$mdyst, coordsts$timest, sep=" ")
coordsts$date_time = ymd_hms(paste(coordsts$mdyst, coordsts$timest))
coordsts$date_hours <- round_date(coordsts$date_time, "hour")
intsph <- coordsts %>%
  select(lat, lon, mdyst, date_hours)
colnames(intsph) <- c("lat", "lon", "date", "hours")
polf <- fortify(pol)
intsph$group <-1

##################### static facet

intsph$date <- as.Date(intsph$date)
intsph$days <- yday(intsph$date)
intsph$days <- intsph$days - 147

p <- ggplot(polf, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour="black", fill=NA, size=.1) + theme_base() +  
  theme_tufte(base_family = "Helvetica") +
  scale_shape_manual("", values=c(8)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank(),
        legend.position = "right",
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 15)) +
  theme(legend.position = "none") + 
  theme(aspect.ratio = 1)
p +
  geom_jitter(data=intsph, aes(x = lon, y = lat), 
              size = 1, alpha = .1, color = "black", width = .2, height = .2) + 
  facet_wrap(~date, ncol=7)

ggsave(file="plots/survfacet.png",
       width = 400, height = 300, units = "mm")